#ifndef __SUPERLU_UTILS_NURBS__
#define __SUPERLU_UTILS_NURBS__

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/storage.hpp>

#include <slu_ddefs.h>

#include "../math/vector.h"

namespace SUPERLU_wrapper {

template<class T, class M>
bool Solve(const boost::numeric::ublas::mapped_matrix<T>& A, const boost::numeric::ublas::vector<T, M>& b,  boost::numeric::ublas::vector<T, M>& x) {
	namespace ublas = boost::numeric::ublas;
	typedef const boost::numeric::ublas::mapped_matrix<T> SparseMatrix;
	typedef boost::numeric::ublas::matrix_column< SparseMatrix > Column;

	std::vector<T> va;
	std::vector<int> vxa;
	std::vector<int> vasub;
	std::vector<T> vb;

	int k = 0;
	for(unsigned j = 0; j < A.size2(); j++) {
		Column c(A, j);
		typename Column::iterator it;
		vxa.push_back(k);
		for(it = c.begin(); it != c.end(); it++) {
			vasub.push_back(it.index());
			va.push_back(static_cast<T>(*it));
			k++;
		}
	}
	vxa.push_back(k);

	for (unsigned j = 0;j < b.size();j++)
		vb.push_back(static_cast<T>(b[j]));	

	int nnz = va.size();
	T* a = va.empty() ? 0 : &va[0];
	int* asub = vasub.empty() ? 0 : &vasub[0];
	int* xa = vxa.empty() ? 0 : &vxa[0];
	T* rhs = vb.empty() ? 0: &vb[0];
	
	// superLU def:	
	SuperMatrix A_, L_, U_, B_;

	int *perm_r;
	int *perm_c;

	int nrhs, info, m, n;

	superlu_options_t options;
	SuperLUStat_t stat;

	// the matrix column and row size
	n = A.size1();
	m = A.size2();

	// create the matrix A
	dCreate_CompCol_Matrix(&A_, m, n, nnz, a, asub, xa, SLU_NC, SLU_D, SLU_GE);

	// Create right-hand side matrix B. 
	nrhs = 1;
	dCreate_Dense_Matrix(&B_, m, nrhs, rhs, m, SLU_DN, SLU_D, SLU_GE);

	// pivoting arrays
	if(!(perm_r = intMalloc(m)) ) {
		std::cerr << "Mallo fails for perm_r[]." << std::endl;
		return false;
	}
	if(!(perm_c = intMalloc(m)) ) {
		std::cerr << "Mallo fails for perm_c[]." << std::endl;
		return false;
	}

	// we can use the strictly diag domiant parameters???
	set_default_options(&options);

	options.ColPerm = NATURAL;

	StatInit(&stat);

	// solve the Ax = B and results will be stored in B
	dgssv(&options, &A_, perm_c, perm_r, &L_, &U_, &B_, &stat, &info);

	// put the value from B to sol:
	DNformat *Bstore;
	register int i;
	T* dp;

	Bstore = (DNformat *) B_.Store;
	dp = (T*) Bstore->nzval;

	for (i = 0; i < B_.nrow; ++i)
		x(i) = dp[i];
	
	// De-allocate storage 
	SUPERLU_FREE (perm_r);
	SUPERLU_FREE (perm_c);
	Destroy_SuperMatrix_Store(&A_);
	Destroy_SuperMatrix_Store(&B_);
	Destroy_SuperNode_Matrix(&L_);
	Destroy_CompCol_Matrix(&U_);
	StatFree(&stat);

	return true;
}

}

#endif //__SUPERLU_UTILS_NURBS__
